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Abstract 


A dynamic model of a 1.2 kW polymer electrolyte membrane (PEM) fuel cell (FC) is developed and validated through a series of experiments. 
This dynamic model is mostly oriented towards control and operation optimization and can be a useful tool for the design of FC-based systems. 
In the methodology proposed, theoretical equations are combined with experimental relations, resulting in a semi-empirical formulation. The 
model assumptions are discussed extensively as the equations are presented. This model contributes to the description of the following areas: fluid 
dynamics in the gas flow fields and gas diffusion layers (oxygen, hydrogen, liquid water and vapor); thermal dynamics and temperature effects; a 
novel algorithm to calculate an empirical polarization curve. As a result, this model can predict both steady and transient states (such as flooding 
and anode purges) due to variable loads, as well as the system start-up. Based on this model, a simulator software package has been developed, 
which is available upon request. The model parameters have been adjusted specifically for a 1.2 kW Ballard stack, which can be considered a 
benchmark as it is widely used by research groups worldwide. Finally, the simulated results are compared to experimental data from the Ballard 


stack, demonstrating the accuracy of the proposed model methodology. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


PEM fuel cells are expected to play an important role in the 
future energy scenario [1], where they will be used in both auto- 
motive and stationary applications, implying that steady-state 
and transient modes will have to be taken into account. The 
importance of having an accurate model that predicts the fuel 
cell behavior is a very important issue. A model not only pro- 
vides a framework for analyzing the performance of the PEM 
fuel cell system and its components but it is also valuable in that 
it can supply the value of internal variables which are difficult 
to measure, such as the water content inside the flow fields. 

There are many PEM fuel cell models in literature. In fact, 
many fuel cell models have been developed over the past 15 
years. Earlier models, such as in Ref. [2], presented an empiri- 
cal polarization curve based on calculated coefficients, as some 
recent papers [3] have shown. In Refs. [4-6] an extended 
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equation, with a larger number of parameters was proposed, 
improving the formulation of the polarization curve depen- 
dency on the stack temperature and the hydrogen and oxygen 
partial pressures. In this paper, an improved equation is pre- 
sented, which properly models the temperature and reactant 
partial pressure influence on the curve, as well as the system 
start up sequence, in conjunction with the dynamic equations 
of the model. Due to the optimized formulation proposed, a 
direct geometric-based algorithm can be deduced, which allows 
a direct equation parameter calculation to be performed using 
four experimental points and two experimental ratios, improving 
the off-line and computationally costly iterative method pro- 
posed in Ref. [4]. This represents an important improvement 
since the equation can now be upgraded on-line and be used, for 
example, in adaptive controllers that could cover the membrane 
degradation as it ages. 

As the polarization curve only includes the steady state by 
itself, later works have focussed on the fluid dynamics inside 
the stack, taking transient behavior into account. Bernardi and 
Vebrunge and Springer et al. [7,8] studied the water flow across 
the membrane and the variable membrane hydration. Gurski 


A.J. del Real et al. / Journal of Power Sources 173 (2007) 310-324 311 


Nomenclature 


species activity 

area (m?) 

mole concentration (mol m~?) 
specific heat capacity (J (kg K)~!) 
heat capacity (J K7!) 

diffusion coefficient (m? s7!) 
effective diffusion coefficient (m? s7!) 
Faraday constant (C mol!) 

mass specific enthalpy (W kg~!) 
mass specific enthalpy of formation (W kg~!) 
enthalpy flow rate (W) 

current (A) 

current density (A m~?) 

valve coefficient (kg (bar s)~!) 
heat transfer coefficient 

mass (kg) 

mass flow rate (kg s7!) 

molar mass (kg mol`!) 
electro-osmotic drag coefficient 
number of fuel cells 

molar flux (mols 7! m7?) 
pressure (N m~?) 

power (W) 

heat flow rate (W) 

ideal gas constant (J (mol K)~!) 
fraction of liquid water volume to the total volume 
level of immobile saturation 
reduced liquid water saturation 
time (s) 

temperature (K) 

volume or voltage (mor V) 

mass fraction 

humidity ratio 

polarization curve coefficient 


a 


a5 28 AA ee a Se n e 


N 
> 


+2 SN Be ROYS Z 


Greek letters 

a) flooding experimental coefficient 

Qw conductivity correction coefficient 

y volumetric condensation coefficient 

ô thickness of diffusion layer (m) 

E porosity or emissivity 

n viscosity (kg (m s)” 1) 

Oe contact angle (°) 

À water content or excess ratio 

u permeability (m7) 

Ur relative permeability (kg mol~!) 

p mass density (kg m~?) 

o surface tension or Stefan—Boltzmann constant 
Nm! or Wm? K~*) 

(0) relative humidity 


Superscripts and subscripts 


a dry air 
act activation 
amb ambient 


anch anode flow channel 
anGDL anode diffusivity gas layer 


atm atmospheric 
B fuel cell stack body 
c capilar 


cach cathode flow channel 
caGDL cathode gas diffusion layer 
conc concentration 

conv convection 

cool coolant 


dry dry 

elec electric 
evap evaporation 
fc fuel cell 
forc forced 

g vapor 

gen generated 


H? hydrogen 
H20 water 


in inlet 

l liquid water 
ma moisture air 
memb membrane 
nat natural 

N2 nitrogen 

(0) initial or reference 
ohm ohmic 

out outlet 

O2 oxygen 

p pore 

purge purge 

rad radiation 
react reaction 

sat saturation 

st fuel cell stack 


[9] also considered the reactant flow control. In Ref. [10], the 
importance of gas hydration was presented, and water trans- 
port and its relation with the membrane thickness, while the 
feeding gases composition was shown in Ref. [11]. Some other 
models have proposed a multi-dimensional study, such as Refs. 
[12-14]. More complex approaches in 3D modeling have also 
been developed [15-19]. Although these contributions are very 
useful for fuel cell design, they require large computational cal- 
culations. Thus, simplified one-dimensional models are more 
suitable for control purposes, such as the ones presented in Refs. 
[20-25], and they allow faster simulations and implementations 
as well. In Ref. [26], both gas diffusion layers (GDLs) and gas 
flow fields are modeled, considering lumped parameters, divid- 
ing each GDL into three control volumes and each flow field 
into one. The work presented here follows that method, but sim- 
plifying the GDLs, considering each one to be a unique control 
volume and thus reducing the computational cost. 
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Concerning thermal dynamics, there are some detailed stud- 
ies such as in Refs. [27—30]. Wetton [31] proposed an explicit 
thermal model to analyze the temperature gradient of different 
layers in the fuel cell stack. Also, Sundaresan [32] presented 
a very detailed 1D thermal dynamic model. In Ref. [33], all 
the principal thermal effects are presented, and then the equa- 
tions are simplified, neglecting the less important coefficients 
and specifically adapting a water cooled stack. Therefore, this 
work is a variation, considering here an autohumidified stack 
with air as coolant. 

Notice that many of the existing models neglect the dynamic 
effects and although some others take these issues into account, 
there is still a lack of validated experiments. The main contribu- 
tion of this work is the development of a dynamic model and its 
validation through real results in a system that is used by many 
research groups. 


Fig. 1. 1.2 kW Nexa power module. 


2. Test equipment 


The experimental data presented in this paper was obtained 
from a 1.2 kW Ballard PEM fuel cell (Nexa Power Module, see 
Fig. 1), which is currently being used by many research groups 
and is representative of the state of the art in PEM technology. 
The benchmark is equipped with a controller which assumes 
the control and safety tasks. The stack is composed of 46 cells, 
each with a 110cm?membrane. The system is autohumidified 
and air-cooled by a small fan. Concerning the hydrogen feeding 
of the fuel cell, a dead-end mode with flushes was adopted. Also, 
a PC was used for the acquisition of the measured values and, in 
order to simulate a variable power demand, the energy produced 
was delivered to an electronic load. 


3. Model description 


In the model methodology proposed, theoretical equations 
are combined with experimental relations, resulting in a semi- 
empirical formulation. The model assumptions are widely 
discussed as the equations are presented. Furthermore, it is com- 
posed of three main modules: electrochemical static model, 
fluid-dynamics model and thermal dynamics model. Charac- 
teristic fuel cell times are shown in Table 1. As can be seen, 
electrochemical dynamics can be neglected as they are several 
orders of magnitude slower than the others, as published in Ref. 


Table 1 
Time constants of the fuel cell dynamics 


Dynamical effect Characteristic time 


Electrochemistry (0) (107!°)s 
Fluid-dynamics 0 (107')s 
Temperature O (107) s 
Content of liquid water O (102) s 


[34]. Finally, the model, though generalized for standard PEM 
fuel cell stacks, contains some parameters which depend on the 
physical dimensions of this system, as well as on several other 
particular issues, such as membrane characteristics. In this paper, 
as presented in Section 2, a 1.2 kW Ballard stack has been specif- 
ically adapted. Thus, the value of all the parameters are presented 
in the subsequent sections as well. 


3.1. Fluid-dynamics equations 


The fluid-dynamics equations consider five control volumes, 
rather than nine as in the case of Ref. [26], therefore reducing 
the number of calculations, while still taking into account all of 
the effects presented in that model. Moreover, additional issues 
have been included, such as water evaporation and condensation 
dynamics. 

The fluid-dynamics block is composed of five interconnected 
sub-blocks, which correspond to: the control volumes of the two 
flow channels; the diffusion gas layers of cathode and anode; the 
transport of chemical species across the membrane. All of these 
control volumes are assumed to be at temperatures equal to that 
of the stack, Ty. The signal criteria adopted depends on the 
direction of the flows, as shown in Fig. 2, where the direction of 
the arrows corresponds to positive values. 


3.1.1. Cathode flow channel 

The inlet air flow is supplied by an air pump and conditioned 
by a humidifier. Thus, the values of the inlet flow sitcach,in, its tem- 
perature Teach,in and its relative humidity @cach,in Must be known. 
Dry air composition will be assumed to be equal to that of the 


ANODE FLOW CATHODE FLOW 
CHANNEL ANODE GDL MEMBRANE CATHODE GDL CHANNEL 
Manch,in Mcach,in 
. m 
M, anGDL2anch v,caGDL2anch 
M42, anGDL2anch 9 caGDL2anch 
Manch,out Mcach,out 


ye 


' 


Fig. 2. Signal criteria. 
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Table 2 

Flow channels parameter set 

Parameters Value 

Seach» anch (m) 1.5 x 1073 
Nfe 46 

AQ. (m?) 110 x 1074 
Veach, Vanch (m?) 7.59 x 1074 
tourge (s) 0.5 
Keachout (kg (bar s)~!) 0.01 
Kanch,out (kg (bar s) 1) 0.001 
Kanch,in (kg (bars) ') 0.2 


atmospheric air. For purpose of simplification, Peach,in = Peach- 
Also, the pressure at the end of the channel equals the atmo- 
spheric pressure, Dcach,out = Patm. The channel is adjacent to the 
cathode gas diffusion layer, with water exchange m], caGDL2cach 
and 7g caGDL2cach and oxygen exchange 710, ,caGDL2cach OCCUT- 
ring between them. Nitrogen exchange will be neglected as it 
is an inert gas. The water flow depends on the difference of 
concentrations between the flow channel and the GDL control 
volumes and is calculated in Section 3.1.3. Due to the signal 
criteria shown in Fig. 2, oxygen flow is negative in the direc- 
tion of the gas diffusion layer. Furthermore, it equals the oxygen 
amount which reacts with the hydrogen. This flow is calculated 
in Section 3.1.3. The differences in oxygen pressure between 
the cathode flow channel and the cathode gas diffusion layer are 
negligible, as proven by works like [35], which show that the 
exchange dynamics due to concentration differences are very 
low. Finally, the values of all the parameters used in the flow 
channel sections are presented in Table 2. 
The inlet air flow is defined as 


1 
™M™Qb,,cach,in = WOd,cach,in M cach, in (1) 
1+ Wr,cach,in 
1 
MN),cach,in = WN),cach,in Mcach,in (2) 
1+ Wr cach,in 
g Wreachin . 
My cachin = 7 ~~~ Meach, in (3) 
1+ Wr,cach,in 
where 
My Pcach,in Psat( Tcach,in) 
Wr,cach,in = (4) 


Ma Pcach,in — Pcach,in Psat(Teach,in) 


considering that the dry air mass fraction equals the 
atmospheric air mass fraction, wo, cach,in = 0.21Mo, M`! 
and WN, cachin = 0.79Mn,Ma', where Ma = 0.21Mo, + 
0.79Mn, ~ 0.02885 kg mol” !. 

Applying mass balance to the control volume and assessing 
the inlet and outlet flows of the channel and the exchange flow 
between it and the gas diffusion layer result in the equations 
shown below. Notice that 7it)cachin = 0, as it can be assumed 
that no liquid water carried by the inlet air enters the cathode 
channel and 7i1},caGDL2cach,in = 0 because the membrane does 
not allow liquid water transport, it only allows gas transport. 
Furthermore, almost all the liquid water condensed inside the 
cathode channel is dragged by the water exhaust, which results 


in dm} cach /dt = 0: 


dm cach g . . F 

de = M\,cach,in — 1,cach,out — Mevap,cach + M\,caGDL2cach 
(5) 

dmy,cach ‘ g ; : 

dt = My cach,in — Mv,cach,out + My,caGDL2cach + Mevap,cach 
(6) 

dio, cach : . . 

ap MOx,cach.in — 10>2,cach,out ~ 109,caGDL2cach (7) 

dmN, cach . A 

a = MN 3,cach,in — MN3,cach,out (8) 

Mma,cach = MO3,cach + MN3,cach + My,cach (9) 


In order to describe the evaporation and condensation dynam- 
ics inside the channel, the equations proposed in Ref. [36] were 
used. The 7iteyap,cach Value is positive when the steam pressure is 
smaller than the saturation pressure, causing water to evaporate 
in that case. When the water condenses, the value is negative. 
Moreover, a logical restriction that would prevent the evapora- 
tion of more water than is available must be considered. In this 
way, 


Mevap,cach 


= min | Age(Psat(Tst) — P v,cach) M1 caGDL2cach 


My 
2a RT 


Pressures inside the channel are calculated as 


Peach = POo,cach + PN3,cach + Pv,cach (11) 
RTs 
POd>,cach = Mo, Voan O2” (12) 
RT 
PNo,cach = MN: Va Oa (13) 
RTs 


Py,cach = Peach Psat(Tst) = My cach (14) 


My Veach 

Supposing that all the liquid water on the surface of the chan- 
nel is dragged by the air that circulates across the cathode, outlet 
flows will be: 


m ma,cach,out = Kcach,out( Pcach — Pcach,out) (15) 
; M™Qp,cach . 
™M™Q),cach,out = ———— ma,cach,out (16) 
M ma,cach 
. MNp,cach . 
MN) ,cach,out = ————™ ma, cach, out (17) 
M ma,cach 
. My cach . 
My cach,out = ————M ma,cach,out (18) 
M ma,cach 
M] cach,out = M] caGDL2cach = Meyap,cach- (19) 
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3.1.2. Anode flow channel 

The equations that model the anode flow channel are anal- 
ogous to the ones that model the cathode channel. In this way, 
applying mass balance, assuming that dry hydrogen enters the 
channel and taking into account the signal criteria shown in 
Fig. 2, result in: 


dim anch . g P 

-u = —M\,anGDL2anch — Mevap,anch — 1,anch,out (20) 
dimy,anch 

ar = —My,anch,out — My,anGDL2anch + Mevyap,anch (21) 
dmH,; anch : R 

— y; tHa.anch.in — FH anch,out — PtH anGDL2anch (22) 


Water steam and hydrogen partial pressure inside the channel 
are defined as 


RT 
Py,anch = M, Vaen VAP (23) 
RTs 
PH2,anch = Mn, Vanor a” (24) 


Condensation dynamics are calculated as in the cathode side: 


Mevap,anch 


= min m 1,anGDL2anch 


(25) 


igota Woa 
fel PsatU st Py,anch ARTs, 


During the experimental time fpurge, most liquid water con- 
densed on the surface of the channel is dragged by the purge gas 
flow, which is calculated as shown below. When the purge valve 
is closed, Vout,open = 0 (see Section 3.4.4), being positive when 
the valve is opened: 


Panch = Py,anch + PH3,anch (26) 

™Mma,anch,out = Kanch,outVout,open(Panch = Patm) (27) 

Mma,anch = MH3,anch + My,anch (28) 

; M¥p,anch . 

MH; ,anch,out = ————— M ma,anch,out (29) 

M ma,anch 

: Myanch . 

My anch,out = —————™ ma, anch, out (30) 
Mma,anch 

: Mlanch .p. 

M],anch,out = if 7M ma,anch,out > 0. (31) 
tpurge 


3.1.3. Cathode and anode gas diffusion layers (GDLs) 

In Ref. [26], each gas diffusion layer is divided into three 
control volumes. The equations presented here simplify that 
study, considering each GDL as a unique control volume. This 
formulation, although being much simpler and less costly com- 
putationally, describes the system behaviour very well. As this is 
a novel approach, almost all the equations included here present 
some departure from those presented in the literature. Partic- 
ularly, the spatial gradients have been linearized, resulting in 
algebraic equations and thus allowing a simpler solution to the 


Table 3 

Gas diffusion layers parameter set 

Parameters Value 

E 0.5 

VopL (m?) 2.53 x 1074 
Dy (m? s7!) 34.5 x 1076 
dcp (m) 00.5 x 1073 
y 0.9 x 103 
Sim 0.1 

o (Nm!) 0.5 

Oe (°) 120 

u m?) 2.55 x 10-133 
m (kg(ms)~!) 978 


problem. Moreover, the parameter values based on those of Ref. 
[35] are shown in Table 3. 

The equations that model all the phenomena occurring inside 
the layers can be divided into two groups: gaseous phase and 
liquid phase. As nitrogen is an inert gas, it will not be taken 
into account. Also, oxygen and hydrogen inside the GDLs will 
be assumed to be at the same pressure as they are in the flow 
channels. Therefore, the hydrogen and oxygen flow between the 
channels and the diffusion layers will be imposed by the electro- 
chemical reaction. Thus, the mathematical formulation is then 
simplified, modeling the water steam concentration gradients, 
the liquid water capillary pressure and the water condensation 
dynamics. 


3.1.3.1. Gas phase. Gaseous species diffusion occurs between 
low concentration regions and higher concentration ones. Hence, 
water steam molar concentrations inside each diffusion gas lay- 
ers are calculated as 


Pv,caGDL 


Cy,caGDL = RT, (32) 
Cyeach = © A (33) 
CyanGDL = a (34) 
ee ae (35) 


The effective diffusion rate, (Dj), is a function of the gas 
diffusion layer porosity [35], € = Vp Vani the fraction of liquid 
water volume to the total volume s, and the diffusion coefficient 
D j: 


V x 
sj= 7 for j= caGDL,anGDL (36) 
P 
e—0.11\ 075 
(Dv,ca) = Dye( =) (1 = Sea)” (37) 
e— 0.11 0.785 
(Dyan) = de( =) (1 = San)? (38) 
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Molar flows are a function of the effective diffusivity and 
concentration gradients: 


Neca a (Dva) (e Saco) (39) 
GDL 
Nyan = (Dan) (2s = arcot | (40) 
ÔGDL 


where ôgpz is the diffusion layer thickness. 
Water steam partial pressures inside the diffusion layers are 
evaluated as 


= = RT (“= = o + Rea ) 
GDL 
(41) 
d Nyan =N 
@Py,anGDL = RT Svan ee 7 Revap,an (42) 
dt ÔGDL 


where Ny memb is calculated in Section 3.1.4 
Evaporation flows are modeled as presented by [35]: 


Psat(Tst) — Pv,caGDL 
R Tst 


(43) 


Revap,ca = 


Psat(Ts) — Pv,anGDL 
RT 


(44) 


Revap,an = 


where, with j = ca,an, a constraint to prevent the evaporation 
of more liquid water than is available must be included: 


ifVj=0 and Revap,j > 0 = Revap,j = 0. (45) 


Lastly, oxygen and hydrogen flows and the amount of water 
steam generated in the electrochemical reaction are calculated 
using stoichiometric balances: 


Ist 


Ne 46 

= on, (46) 
Tet 

No» react = TFA; (47) 
Is 

Nu react = Ir (48) 


F being the Faraday constant. 
Mass gaseous flows exchanged between gas diffusion layers 
and flow channels are calculated below: 


MH, anGDL2anch = Afcnfe Mn, NB ,react (49) 
MO,,caGDL2cach = Afcnfe Mo, No,,react (50) 
My anGDL2anch = Åfcnfe My Nyan (51) 
My caGDL2cach = AfeNfc My Ny ca- (52) 


3.1.3.2. Liquid phase. Liquid water volume evaluation is based 
on mass balances in both the anode and the cathode sides, as 


dVi caGDL 
A ~ 


I ~ —M],caGDL2cach — Revap,ca Mye VGDL (53) 


Table 4 
Membrane parameter set 
Parameters Value 
Qw 15 
memb (m) 35 x 1076 
Pmemb, dry (kg m7?) 2x 103 
M memb,dry (kg mol” !) 1.1 
dVianGDL 
ay ManGDL2anch — Revap,anMyé VGDL (54) 


where p is the liquid water density. 
Let j = ca,an. The reduced liquid water saturation S; j, Sim 
being the liquid water immobile saturation [35], is modeled as 


>j im ara <s;<1 
Si Sj 
S27) (ag, a (55) 
0 paraO < s; < Sim 


Furthermore, the capillary pressure pe is calculated via the 
Leverette J function, which describes the relationship between 
the capillary pressure and the reduced liquid water saturation S,. 


o COS be 


Pe = Kje” 


[1.417S; — 2.1208? + 1.26352] (56) 
——qy—_q— 
J(Sr) 


where o is the surface tension which corresponds to water and 
air, 0, is the contact angle and u is the absolute permeability. 

At this point, the capillary flows of liquid water rn caGDL2cach 
and M1 anGDL2anch Can then be defined as in Refs. [35,37], where 
Urn = sS is the relative permeability of liquid water and 7 is 
its viscosity: 


: AfcnfcHHr | dpe] Sca 
M\,caGDL2cach = (57) 
nı dS | égpi 
; AfcnfeHHr |dpe| S: 
M1 anGDL2anch = —_ =| (58) 
nı dS | égpi 


3.1.4. Membrane 

As the equations presented in this section are analogous to 
the ones published in Ref. [26], the main innovation regarding 
the membrane is the parameter adaptation for the Ballard stack, 
which is presented in Table 4. 

The membrane, being waterproof, does not allow the circu- 
lation of liquid water but does permit gas diffusion. Therefore, 
oxygen, hydrogen, nitrogen and water steam are exchanged 
between both sides of the membrane. As the anode is fed with 
pure hydrogen, nitrogen flow occurs in the direction of the anode. 
This phenomenon is negligible due to its very slow dynamics 
and also because nitrogen is an inert gas that does not affect 
the electrochemical reaction. Moreover, the amount of nitrogen 
inside the anode will never be significant, as it is dragged during 
purging. Water steam flow is a more complex issue as it is a 
convergence of two distinct phenomena: ‘electro-osmotic drag’ 
and ‘back diffusion’. Thus, the molar flow across the membrane, 
published by Springer et al. [38], is described as 


Cy,ca — Cv,an 


(59) 


dw Dy 


Ny memb = Ma 
F memb 
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where j = Ist/Afe is the current density, ng the electro-osmotic 
drag coefficient, Dy the mass diffusivity of water vapor in the 
membrane, dmemb the membrane thickness and cy cq and Cyan are 
the water concentrations in both sides of the membrane. œw is an 
experimental parameter which corrects the possible deviations 
of the experimental values obtained from the literature, as they 
might be obsolete due to their being based on older membranes. 

The concentrations cy, j, where the subscript j corresponds to 
both the cathode and the anode, are calculated as 


Pmemb,dry A (60) 


Cy j = 
j Mmemb,dry 

where pPmemb,dry 18 the dry membrane density, Mmemb,dry the dry 
membrane molecular weight and A; is calculated as 


Aj = 0.043 + 17.81aj; — 39.854; + 36.04; (61) 


aj being equal to ¢jGp_, the relative humidity of the gas inside 
the gas diffusion layers. 

Dy is calculated as the piece-wise linear approximation 
shown below: 


1 1 
Dy = D 2416 | — — — 62 
+= Pian (2416 (555-72) “= 

1g", Aan < 2 
107101 + 2(Aan — 2)), 2<A <3 

D, = a (Aan — 2)) <i< (63) 
107103 — 1.67(Aan — 3), 3 <A < 4.5 
1.25 x 10719, han > 4.5 


where j corresponds to an,ca and D,,, is the corrected diffu- 
sivity coefficient. Lastly, the electro-osmotic drag coefficient, 
described by Dutta et al. [15], is calculated using: 


na = 0.002922, + 0.05ran — 3.4 x 107). (64) 
3.2. Electrochemical equations 


This section presents one of the main contributions of this 
paper. In fact, in the mathematical description of the polarization 
curve, despite being based on the same general equation as the 
one extensively used in the literature (see Eq. (65)), all of the 
equation terms are presented in a novel way, thereby optimizing 
and simplifying its formulation. As a result, a direct algebraic 
algorithm can be geometrically deduced, in order to adapt the 
polarization curve so that it corresponds to the experimental data 
associated with a given PEM stack. All of these parameters, 
which are obtained below in this section for the Ballard stack, 
are shown in Table 5: 


Vie = VO — Vact — Vohm — Veonc- (65) 


The voltage supplied is evaluated by curves that present the 
voltage of the cell versus the current density, j = I,,/ Age. Fig. 3 
shows different real voltage data corresponding to several values 
of temperature, oxygen partial pressure and current density taken 
from the Ballard stack, while the hydrogen partial pressure was 
maintained nearly constant and equal to 1.25 bar. Notice that 
the upper graph shows voltage data for different current density 


Table 5 

Electrochemical parameter set 

Parameters Value 

(Pii, Piv) (0, 1) 

(Pai, pay) (0.06, 0.785) 

(P3i, P3v) (0.4, 0.555) 

(Pai, Pav) (0.5, 0.35) 

Paves 0.16 (bar) 

Plian 1.25 (bar) 

T? 308 (K) 

AVic —3 -l1 

ate 2.93 x 10-3 (VK~') 
ANg 7.61 x 107! (V bar™!) 


APO3,ca 
(x1, X2, X3, X4, X5, X6, X7, X8) (1.17,7.61 x 1073, 0.24, 0.18, 1.50 x 
10-7, 0.64, 288.59, 10.00) 


æ 15 


values, while the middle graph presents the temperature of the 
stack for each experimental point. Lastly, the lower graph shows 
the calculated oxygen partial pressure corresponding to each 
value. The main factors that influence the polarization curves 
are: 


e Current: open circuit voltage (vo) falls as the current sup- 
plied by the stack increases. Thus, in the first stage, up to a 
certain current value, activation overvoltage drops (Vact) pre- 
vail, as aresult of the need to move electrons and to break and 
form chemical bonds. At a later stage, as current density rises, 
ohmic losses (Uphm) prevail. They are derived from membrane 
resistance to transfer protons and from electrical resistance of 
the electrodes to transfer electrons. When current is very high, 
at maximum power level, concentration overvoltage (Veonc) 
produces a quick drop of the voltage due to internal ineffi- 
ciencies at high levels of reactives consumption. Although 
these last drops are modeled by the equation proposed here, 
such high current values are not normal as they can cause fast 
degradation of the membranes of the cells. In Fig. 4, the polar- 
ization curve that corresponds to the 1.2 kW stack modeled 
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Fig. 3. Experimental data for several partial pressures and temperatures. 
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Fig. 4. Fuel cell polarization curve and voltage drops contributions. 


here can be seen. The drops described above are case-adapted 
as well. 

Reactant partial pressures variances: oxygen partial pres- 
sure inside the cathode and hydrogen partial pressure inside 
the anode notably influence the voltage supplied by the fuel 
cells. As the pressure of any of these gases rises, the voltage 
increases for every current level. Due to the typical dead- 
end configuration of fuel cell stacks system, hydrogen partial 
pressure inside the anode remains nearly constant. Thus, the 
oxygen partial pressure is that which fluctuates more during 
power demand changes, as it depends on auxiliary equipment 
dynamics. The oxygen influence can be seen in Fig. 5. The 
surface shown was obtained solving the equations proposed 
and plotting the polarization curves corresponding to values 
of oxygen partial pressure ranging from 0 to 0.25 bar, while 
maintaining the stack temperature and the hydrogen partial 
pressure constant at 310K and 1.25 bar, respectively. Also, 
experimental data was included, thus comparing theoretical 
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Fig. 5. Oxygen partial pressure influence on the calculated polarization curve 
and experimental data (x). 
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Fig. 6. Temperature influence on the calculated polarization curve and experi- 
mental data (x). 


and measured values. Notice the logarithmical tendency of 
the curves, which is described below in this section, partic- 
ularly at low oxygen partial pressures, allowing the fuel cell 
system start-up sequence to be modeled (see also Fig. 8). 
Temperature: fuel cell stack temperature has many effects on 
fuel cell performance, including changing the activity of the 
catalyst, the membrane humidification state, the saturation in 
GDLs and gas diffusion. For purposes of simplification, we 
will assume that, on one hand, temperature changes affect 
the gas pressure inside both anode and cathode. On the other 
hand, membrane characteristics change as the temperature 
rises, which results in a voltage increase at every current level. 
As a result, it is important to distinguish between these two 
effects when realizing an experimental characterization of the 
fuel cells. Due to the technical constraint of including pressure 
sensors inside the flow channels, mathematical modeling is 
important to obtain simulated values of these pressures, as 
well as to discriminate the two effects described above. Fig. 6 
shows the calculated polarization curve and its dependency 
on temperature (ranging from 280 to 340 K), maintaining the 
oxygen partial pressure at a constant value of 0.15 bar and the 
hydrogen partial pressure at a value of 1.25 bar. Experimental 
data was included to show the accuracy of the theoretical 
model proposed here. 

Amount of condensed water: liquid water occupies pore vol- 
ume and reduces the gas diffusion layer surface resulting in 
voltage and fuel cell efficiency drops. This phenomenon is not 
important in the cathode side as the water excess is dragged 
by the air flow, but water accumulates in the anode side until 
it is dragged during a purge event. This effect has been mod- 
eled around the assumption that as the amount of liquid water 
inside the diffusion gas layer increases, the membrane effec- 
tive area decreases, resulting in an increased current density 
and hence a drop in voltage, as shown in Fig. 7, where the cal- 
culated current density j = Ist/Afce effective from the measured 
value of Zst and the calculated value of A fc effective is presented 
in the upper graph, while the measured and calculated voltages 
are plotted in the lower graph. Some papers, such as in Ref. 
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Fig. 7. Current density during purge events. 


[26], present also the same approach, as under testing con- 
ditions, they have proven experimentally that cathode purges 
have little effect on cell voltages. At the same time, they found 
that voltage significantly recovered following an anode purge. 
This phenomenon can be explained since the liquid water 
condensed in the cathode side is continuously dragged by the 
exhaust air crossing the cathode. As the anode is in dead-end 
mode, the liquid water accumulated is not dragged, except 
after purges. Also, as discussed in Ref. [26], this water forms a 
thin film (experimentally measured in Ref. [39]) blocking part 
of the active fuel cell area (resulting in a lower apparent active 
fuel cell area value A fc effective) and thus increasing the current 
density. 


Thus, the following equation is proposed to model all the 
aspects described above: 


Vic = x1 + 2(Tot — T2) + x3(0.5 In(po,,ca) + In(pu,)) 
KS ee eS 


vo AVye/ATot AVic/Ap 
xa(1 — exp(—j/xs)) — x6» j — x7 < jt) (66) 
N ———— ier Sa 

Vact Vohm Vconc 


D 7 
Afc Ag. — QM ),anch) 
Se 


(67) 


Afe effective 


The logical constraint of Vi. being positive or zero, must be 
included with the stack voltage, being defined as the sum of all 
individual cell voltages: 


Vic = 0 (68) 


Vst = Nfe Vee. (69) 
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Fig. 8. Oxygen pressure in different operation ranges. 


The assumptions made in order to obtain the electrical equa- 
tion are: 


e Distinction between temperature and oxygen and hydrogen 
pressure influences. 

e Linear voltage variations near a nominal temperature en 
This does not imply less generality as it can be experimen- 
tally proven that linear range covers almost all the normal 
operational range of the stack. 

e Logarithmic voltage variations as reactant partial pressures 
change. This allows the start-up sequence modeling. How- 
ever, notice that linearity can be assumed near the nominal 
points Ponca and Pan: On one hand, this assumption makes 
the determination of A Vfc/ Apo, ,ca much easier. On the other 
hand, it simplifies the algebraic algorithm used to obtain 
the parameters. The logarithmic tendency is more important 
during start up and shut down stages, when reactant partial 
pressures are outside of the nominal operational range (see 
Fig. 8). 


The purpose of the algebraic algorithm presented below is 
the calculation of the parameters x1, ..., xg. Four experimen- 
tal points, (pji, pjv) for j = 1,...,4, which define the shape 
of the polarization curve, are necessary. Also, the coefficients 
AVic/ATs and A Vfc/Apo,,ca Must be known. The four points 
must be chosen as they divide the polarization curve into three 
parts, each one corresponding to one of the three voltage drops 
described at the beginning of this section. Also, A Vice /ATgt and 
A Vfc/APo;,ca are experimentally obtained from the operation 
of the fuel cell stack near the nominal points T? and Poses 
respectively. In this way, and after scaling the experimental data 
shown in Fig. 3 using the values AV. /ATy and AV¢e/ Apo, ca, 
a polarization curve at a fixed temperature (308 K) and oxygen 
(0.16 bar) and hydrogen (1.250 bar) partial pressures is obtained. 
In Fig. 9 an example of a possible choice of the four points 
needed by the algorithm is shown. Also, the experimental data 
were collected after a purge event, so the measurements were not 
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Fig. 9. Point choice for the parameter identification algorithm of the polarization 
curve. 


affected by the flooding effect. Comparing Figs. 3 and 9, data 
alignment can be observed. Notice that this can be taken as an 
experimental verification of all the assumptions of linearity and 
distinction between temperature and reactant partial pressures 
discussed above: 


1+ pt. 
teks t+ Päi (10) 
0.25 paj 
_ (Pav — Pav) + (Pow — P3v)((Pai — P3i)/(P3i — P2i)) 
— phi? + p$t*®) (pai — psi)/(p3i — p2i)) 
(71) 
= = (1+) 
ns (Poy — P3v) — x7 P3i (72) 
(p3i — P2i) 
maS P2i z Pli (13) 
X4 = Piv — Pov — X6 P2i (74) 
AVi 
0) Cc 
= 75 
X3 Po, ca Ap T ( ) 
A Vfe 
= 76 
X2 AT, (76) 
x1 = Py — X3(0.5 (pO, ca) + In(Pe, an))- (77) 


Lastly, a is obtained measuring the experimental voltage 
drop that occurs as liquid water accumulates in the anode. Due 
to the technical constraint of placing a sensor inside the anode 
to measure the amount of liquid water accumulated there, an 
estimated value is then necessary, resulting ina} = 15, as shown 
in Table 5. 


3.3. Thermal equations 


As published in Ref. [33], a simplified thermal model is pro- 
posed here, taking into account only the main terms of the overall 


Table 6 

Thermal parameter set 

Parameters Value 
mg (kg) 5 

Cs kg! KT!) 1100 
E 0.9 
Ap2Amb,ad (m°) 0.1410 
AB2Amb,conv (m?) 0.0720 
AB2Amb,conv (m?) 1.2696 
hB2Amb,nat (W kg™!) 14 

Kri 0.0156 
Ky 1 


energy balance. This assumption results in an easier experimen- 
tal adjusting of the equation parameters. The contribution of this 
paper to this field is the adaptation of the equations proposed in 
Ref. [33], so that then characterize the 1.2 kW Ballard stack. As 
discussed before, the stack modeled in Ref. [33], as in almost 
all other literature, corresponds to a water-cooled stack system. 
Here, however, an air-cooled stack has been modeled, and the 
physical parameters have been obtained as well (see Table 6). 

An energy balance is done in order to obtain the thermal 
model, taking into account the energy produced in the chemical 
reaction of water formation (which is supposed to be formed as 
water steam) Hreac, the energy supplied in the form of elec- 
tricity Pelee and the amount of heat evacuated by radiation 
Orad,B2amb and both natural and forced convection Ü conv,B2amb. 
Heat removal is completed through forced convection by a small 
fan. In bigger fuel cell stack systems, where the amount of heat 
is considerably larger, water cooling is necessary. In those cases, 
the forced convection term should be substituted by other terms 
which model heat exchange in cooling fluid. The energy balance 
results in: 

dT 

MstCst— = Hyeac — 


dt P, elect — Ü rad.B2amb (78) 


The enthalpy flow rate, where hv 0@) is the mass spe- 
cific enthalpy of formation of water steam and cp,H,, Cp,O, and 
Cp,H20(g) are the specific heats of hydrogen, oxygen and water 
steam respectively, will be: 


Hreac = MH reac Alp, + MoO, reac Aho, 


= Wu0,g0n(2) (AG ol) Ahm oe) (80) 
Ahy, = Cp, (Tanch,in — T?) (81) 
Aho, = Cp,03(Teach,in — T’) (82) 
Ahy, ole) = Cp, H20(e)(Tst — T°) (83) 


where 7° is the reference temperature for the enthalpy. 
Energy yielded in the form of electricity is calculated as 


Palec = Vst Ist (84) 


Heat exchanged as radiation, where ¢ is the emissivity, o = 
5.678 x 10-8 W m~? K~‘, the Stefan-Boltmann constant and 
AB2amb,rad İS the radiation exchange area, is modelled as shown 
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Fig. 10. Nexa power module system. 


below: 


OradB2amb = ECAB2amb,rad(T¢ — T$a) (85) 


Lastly, the convective term is composed by two others, one 
of them corresponding to the natural convection O conv, B2amb,nat 
and the other, to the forced convection O conv,B2amb, forc- In each 
case, the convective heat transfer coefficients (AB2amb,nat and 
hp2amb,forc) are different, just as the exchange areas (AB2amb,conv 
and Ap?cool,conv) are, because natural convection takes place in 
the fuel cell lateral walls, and forced convection occurs across 
the internal walls of the cells, which are constructed as a radiator. 


Oconv,B2amb = Qconv,B2amb,nat +F Q conv,B2amb, forc (86) 
Q conv,B2amb,nat = (AB2amb,nat AB2amb,conv )(Tst = Tamb) (87) 
Q conv,B2amb, forc = (hB2amb,forc AB2cool,conv (7st = Teool) (88) 


where 


hp2amb,forc = Kpi(Weoot)*®?. (89) 


3.4. Auxiliary equipment 


Fuel stacks are formed by numerous fuel cells connected in 
series, but in order to supply energy, they need auxiliary equip- 
ment. These devices differ according to the application they are 
designed for and fuel cell stack size. In this work, a commercial 
1.2kW Nexa module power system is modeled. In Fig. 10, a 
diagram of the auxiliary equipment can be seen. The devices 
that are particularized are: air pump, cooling fan, humidifier and 
inlet hydrogen conditioning system. 


3.4.1. Air pump 

The air pump is modeled as a black box. Its transfer function, 
which is identified using experimental data obtained during a 
test, is presented below. A transfer function is a mathematical 
representation of the relation between the input and output of a 
system, and is commonly used in signal processing, communi- 
cation theory, and control theory. A model with the following 
discrete time linear transfer function (or transfer operator, see 


e.g. [40]) is used: 


0.06229z7 — 0.0145z 
u 
z? — 2.27722 + 1.811z — 0.5016 


where z denotes the shift operator, i.e. y(t) = y(t — 1)z and t 
stands for the current discrete time instant. Here, the model input 
u(t) is the air pump voltage, which is an analog value ranging 
from 0% to 100%, corresponding 100% to the maximum air 
pump power, and the model output y(t) is the air mass flow 
supplied to the fuel cell stack (kg s~!). In Fig. 11 a comparison 
between the modeled and the real system can be seen. 


y(t) = (t) — 45 (90) 


3.4.2. Humidifier 

Humidifier design is a very complex issue about which there 
is no available information, as their designers are very confiden- 
tial regarding new developments. In the model proposed here, the 
water formed due to the chemical reaction inside the fuel cells is 
utilized to humidify the inlet air flow. Although some mathemat- 
ical approaches can be applied, there is insufficient experimental 
data to validate these models. Thus, it will be assumed that the 
inlet air flow is optimally humidified after passing across the 
humidifier, and there is always enough liquid water available 
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Fig. 11. Air pump real and simulated behavior. 
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from the chemical reaction. This approximation can be done as 
it is very close to the real operational conditions. 


3.4.3. Cooling fan 

In this application, cooling is performed by forced convec- 
tion. A small fan is used to supply the cooling air flow. The 
thermal dynamics of the fuel cell are several magnitude orders 
lower than the fluid-dynamics associated with the cooling air 
flow, meaning that these dynamics are negligible. Moreover, the 
amount of air supplied by the fan can be considered as linearly 
proportional to the control signal of the fan. In this way, the 
equation that links the model input u(t), which is the fan volt- 
age, between 0% and 100%, with the model output y(t), which 
is the air flow supplied expressed in kg s~!, is calculated as 


y(t) = 36u(t). (91) 


3.4.4, Hydrogen feed system 

As can be seen in Fig. 10, the fuel cell system is composed 
of various pressure regulators and valves that enable or disable 
the hydrogen supply. The inlet manifold can be modeled as a 
regulatory valve that adapts the high hydrogen pressure of the 
tanks to a more suitable and lower pressure inside the anode. 
Due to this valve, the inlet hydrogen flow equals the hydrogen 
consumed in the chemical reaction. Thus, vin open being equal 
to 0 when the valve is closed and | when it is opened, Panch 
being the pressure inside the anode, Ptank being the tanks pres- 
sure, 7H ,an,in being the hydrogen inlet flow and Kan,in being a 
characteristic constant of the valve, the model will be: 


MH; ,an,in = Kanch,in Vin,open(Panch — Prank) (92) 
In the same way, the purge valve is modeled as 
Mma,an,out = Kanch,outVout,open(Patm)- (93) 
4. Results and discussion 
After exposing the model equations, implementing them in 


Matlab/Simulink and particularizing the physical and experi- 
mental parameters of the model to the 1.2kW Nexa power 


module, the last step is the model validation. The differential 
equations have been solved using the ode23s solver method, 
provided by the cited Matlab/Simulink software package. As 
have been discussed, the fuel cell stack and its auxiliary devices 
(which are: air pump, humidifier, cooling fan and hydrogen feed 
system) have been modeled. Therefore, the inputs to the model 
are the air pump voltage, the fan voltage, the ambient temper- 
ature, the anode purge valve status and the current demanded 
from the fuel cell system (see Fig. 12). Notice that, depending 
on the current demand, the onboard controller of the fuel cell 
system adjusts the air pump and the fan voltages, as well as 
the hydrogen purge value (which opens to recover the voltage 
which falls due to the flooding effects). Thus, real data from 
tests, collected during several experiments, would be the model 
inputs. The system outputs, which are linked with the inputs 
through the model equations, are the fuel cell stack voltage and 
its temperature, and are also measured. The validation process 
consists of simulating the model applying the real system inputs 
as model inputs, and then comparing the system outputs with 
the simulated ones. In this way, four experiments are proposed: 
start-up stage, constant load (to observe the flooding effect), 
variable load (to analyze the transitory effects) and a long dura- 
tion test with load changes (to validate the thermal equations of 
the model). 


4.1. Start-up sequence 


During the ‘stand by’ stage, both the air pump and the cooling 
fan are off and the hydrogen inlet valve is closed. The period 
required for nominal conditions to be reached (which happens 
when the voltage reaches a nominal stationary value), is known 
as the start-up sequence. When starting up, the air pump and the 
cooling fan switch on and the hydrogen inlet valve opens. At 
the same time, the purge valve stays opened (for approximately 
15 s) in order to ventilate the anode channel. Notice that, during 
this time, no current is delivered by the fuel stack, so there is 
no hydrogen consumption. When a value of approximately 40 V 
is reached, the start up sequence finishes, and power can then 
be delivered. As there is no power demand, anode flooding is 
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Fig. 13. Input data during the start up sequence. 


not serious and therefore cannot be detected when the purge 
valve opens. Real input and output evolution can be seen in 
Figs. 13 and 14, as well as the simulated voltage and temperature 
outputs. 


4.2. Constant load: flooding effects 


During normal fuel cell operation, liquid water condenses and 
accumulates in the anode as it is in dead-end mode, decreasing 
the efficiency and thus requiring anode purges. Thus, the anode 
outlet valve purges for approximately 2s. At the same time, as 
discussed before in this paper, the liquid water accumulated in 
the cathode is removed continuously by the cathode air exhaust, 
so there is no need for cathode purges. The experimental data 
presented was collected applying a constant load, showing the 
flooding effect as well as the voltage recovery after an anode 
purge (see Figs. 15 and 16). Notice that when the load changes, as 
in the other experiments presented in this section, anode flooding 
occurs as well, but as the voltage drop caused by the flooding is 
not very high and voltage is continuously changing due to current 
demand variations, it is difficult to observe flooding effects in 
those cases. This is the reason why a constant load experiment 
was performed in order to show the voltage recovery after a 
purge event. 
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Fig. 14. Output data during the start up sequence. 
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Fig. 16. Output data during flooding test. 
4.3. Variable load: transitory effects 


During abrupt load changes, the current delivered by the fuel 
cell stack varies. At that point, the fuel cell system onboard con- 
troller adjusts the air pump and the cooling fan voltages, in order 
to maintain appropriate air feeding and stack temperature, as 
well as the purge valve, in order to mitigate the flooding effects. 
These data, including also the ambient temperature, can be seen 
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Fig. 17. Input data during transitory effects test. 
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Fig. 18. Output data during transitory effects test. 


in Fig. 17, which corresponds to the model inputs. Solving the 
model equations, the simulated values of voltage and tempera- 
ture can then be compared with the measured ones, as seen in 
Fig. 18. Notice that, when the power demand increases, the air 
supply equipment takes some time to provide the new amount of 
air required, as the electromechanical dynamics present various 
orders of magnitude higher than fluid dynamics response time. 
During this period, the oxygen partial pressure tends to decrease 
and thus the voltage of the fuel cell stack decreases. Also, during 
the test, an anode purge can be observed, but notice that, as it 
was performed during a load change, it is difficult to detect the 
voltage recovery. 


4.4. Variable load: temperature 


As thermal effects are much slower than others, a longer test 
is needed to validate the set of equations dedicated to the stack 
temperature. In the same fashion as the previous experiments, 
load changes were applied to the fuel cell system, collecting 
then the sensor measurements and afterwards, solving the model 
equations. As can be seen in Figs. 19 and 20, the simulated 
temperature tracks the real data thoroughly. Also notice that, 
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Fig. 19. Input data during thermal effects test. 
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Fig. 20. Output data during thermal effects test. 


even when many anode purges are performed, as the graph scale 
is very large and the voltage recoveries are small, it is difficult 
to observe them. 


5. Conclusions 


This work has presented a complete dynamic model of a fuel 
cell. It is a control-oriented model, that can be used for con- 
troller design and optimal operational strategies development 
for FC-based power systems. The model has been validated on a 
Ballard 1.2 kW PEM fuel cell which can be considered a bench- 
mark, since it is well-known by many research groups as a good 
example of the state of the art in this technology. The proposed 
model’s main contributions are related to a novel way of exper- 
imentally obtaining the polarization curve, to the adaptation of 
the thermal equations for an air cooled stack and the modeling of 
the flooding event. Moreover, the model has been validated on a 
real plant, showing that simulated data match experimental data 
quite closely. The simulator software is available upon request. 
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